//
// GPSMapEdit
// (c) Konstantin Galichsky (kg@geopainting.com), 2007
//
// Mercator projection.
//

# include "StdAfx.h"
# include "Mercator.h"
# include "Ellipsoid.h"
# include "PtAlgo.h"
# include "Globals.h"

CMercatorToLL::CMercatorToLL (double _a, double _e, double _k0, double _FE, double _FN, double _x0) :
	m_a (_a),
	m_e (_e),
	m_k0 (_k0),
	m_FE (_FE),
	m_FN (_FN),
	m_x0 (_x0)
{}

// See http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34f.html
void CMercatorToLL::ApplyTo (point_t & _point) const {
	_point.x = (_point.x - m_FE)/(m_a*m_k0)*180/c_PI + m_x0;

	const double t = ::exp ((m_FN - _point.y)/(m_a*m_k0));
	const double chi = c_PI/2 - 2*::atan (t);
	const double e2 = m_e*m_e;
	const double e4 = e2*e2, e6 = e2*e4, e8 = e4*e4;
	const double phi = chi + (e2/2 + 5*e4/24 + e6/12 + 13*e8/360)*::sin (2*chi) +
		(7*e4/48 + 29*e6/240 + 811*e8/11520)*::sin (4*chi) +
		(7*e6/120 + 81*e8/1120)*::sin (6*chi) + (4279*e8/161280)*::sin (8*chi);

	_point.y = phi*(180/c_PI);
}

CLLToMercator::CLLToMercator (double _a, double _e, double _k0, double _FE, double _FN, double _x0) :
	m_a (_a),
	m_e (_e),
	m_k0 (_k0),
	m_FE (_FE),
	m_FN (_FN),
	m_x0 (_x0)
{}

void CLLToMercator::ApplyTo (point_t & _point) const {
	const double lambda = _point.x*c_PI/180;
	const double phi    = _point.y*c_PI/180;
	const double fSinPhi = ::sin (phi);

	_point.x = m_FE + m_a*m_k0*(lambda - m_x0*c_PI/180);
	_point.y = m_FN + m_a*m_k0*::log (::tan (c_PI/4 + phi/2)*::pow ((1 - m_e*fSinPhi)/(1 + m_e*fSinPhi), m_e/2));
}

# ifdef _DEBUG
void Test_Mercator () {
	// Makassar / NEIEZ
	CEllipsoid eBessel1841 ("", 6377397.155, 1/299.15281);
	CMercatorToLL m2ll (eBessel1841.a, ::sqrt (eBessel1841.e2), 0.997, 3900000.00, 900000.00, 110);
	CLLToMercator ll2m (eBessel1841.a, ::sqrt (eBessel1841.e2), 0.997, 3900000.00, 900000.00, 110);

	point_t p (5009726.58f, 569150.82f);
	m2ll.ApplyTo (p);

	assert (::fabs (p.x - 120) < 1e-5);
	assert (::fabs (p.y + 3  ) < 1e-5);

	ll2m.ApplyTo (p);

	assert (::fabs (p.x - 5009726.58f) < 1e-5);
	assert (::fabs (p.y - 569150.82f)  < 1);
}
# endif

